"""
@Author: ISSCA_NEW
@Time: 2024/12/7 16:12
@User: DELL
@File: match.py
@Software: PyCharm
@Attention: 安耐毁誉，八风不动
            利、衰、毁、誉、称、讥、苦、乐
"""
from Tools import base_tools as bt
from Tools import mathtools as mt
# from Tools import gnss_tools
import numpy as np
import os
import netCDF4 as nc
import datetime
import math
import time


class write_nc:

    def __init__(self, filename: str):
        self.ncfile = nc.Dataset(filename + ".nc", "w")

    def write_one_list(self, x: list, name: str, data_type: type or str):
        x_dim: nc.Dimension = self.ncfile.createDimension(name, len(x))
        x_var: nc.Variable = self.ncfile.createVariable(name, data_type, (name))
        x_var[:] = x

    # def write_ddm_value(self, x: list or np.array, name: str, data_type: type or str):
    #     counts = self.ncfile.createDimension(name + "_" + "counts", len(x))
    #     delay_raw = self.ncfile.createDimension(name + "_" + "delay", 17)
    #     doppler_raw = self.ncfile.createDimension(name + "_" + "doppler", 11)
    #
    #     counts_var = self.ncfile.createVariable(name + "_" + "counts", data_type, (name + "_" + "counts",))
    #
    #     delay_var = self.ncfile.createVariable(name + "_" + "delay", data_type, (name + "_" + "delay",))
    #     doppler_var = self.ncfile.createVariable(name + "_" + "doppler", data_type, (name + "_" + "doppler",))
    #     data_var = self.ncfile.createVariable(name, data_type, (
    #         name + "_" + "counts", name + "_" + "delay", name + "_" + "doppler"))
    #
    #     data_var[:] = x

    def close(self):
        self.ncfile.close()

def az(x1:float,y1:float,x2:float,y2:float)->float:
    angle = 0
    dy = y2 - y1
    dx = x2 - x1
    if dx == 0 and dy > 0:
        angle = 0
    if dx == 0 and dy < 0:
        angle = 180
    if dy == 0 and dx > 0:
        angle = 90
    if dy == 0 and dx < 0:
        angle = 270
    if dx > 0 and dy > 0:
        angle = math.atan(dx / dy) * 180 / math.pi
    elif dx < 0 and dy > 0:
        angle = 360 + math.atan(dx / dy) * 180 / math.pi
    elif dx < 0 and dy < 0:
        angle = 180 + math.atan(dx / dy) * 180 / math.pi
    elif dx > 0 and dy < 0:
        angle = 180 + math.atan(dx / dy) * 180 / math.pi
    return angle

class matched:
    def __init__(self):
        self.content:dict={
            "lat":[],"lon":[],"ews":[],"ewsu":[],"ewsv":[],"iws":[],"sst":[],"skt":[],"swh":[],"nbrcs":[],"les":[],"snr":[],
            "inc":[],"cynum":[],"channel":[],"scvx":[],"scvy":[],"scvz":[],"txvx":[],"txvy":[],"txvz":[],
            "azi":[],"ele":[],"prn":[],"sp_az_body":[],"sp_az_orbit":[],"sp_rx_gain":[],"a1":[],"a2":[],"edgea":[],
            "ews100":[],"ews100u":[],"ews100v":[],"eirp":[],"time":[],"sr":[],"pr":[],"nmi":[],
        }
    def save(self,out_file):
        out_nc=write_nc(out_file)
        for parameters in self.content.keys():
            out_nc.write_one_list(self.content[parameters],parameters,"f")
        out_nc.close()




def convert_lon(lon:float)->float:
    if lon>=0:
        return lon
    else:
        return 360+lon

def spherical_coordinates(x: float, y: float, z: float) -> list:
    azimuth: float = np.arctan2(y, x)
    azimuth_deg: float = np.degrees(azimuth)
    azimuth_deg: float = azimuth_deg if azimuth_deg > 0 else 360 + azimuth_deg
    horizontal_dist: float = np.sqrt(x ** 2 + y ** 2)
    elevation: float = np.arctan2(z, horizontal_dist)
    elevation_deg: float = np.degrees(elevation)
    return [azimuth_deg, elevation_deg]

def angle_between_vectors( vec1: np.array or list, vec2: np.array or list) -> float:
    vec1: np.array = np.array(vec1)
    vec2: np.array = np.array(vec2)
    dot_product: np.array = np.dot(vec1, vec2)
    norm_vec1: np.array = np.linalg.norm(vec1)
    norm_vec2: np.array = np.linalg.norm(vec2)
    cosine_similarity: np.array = dot_product / (norm_vec1 * norm_vec2)
    angle_rad: float = np.arccos(np.clip(cosine_similarity, -1.0, 1.0))
    return angle_rad

def topocentric_coordinate_system(sp_loc: list or np.array, rec_loc: list or np.array, sp_lon: float,
                                  rec_lon: float) -> list:
    A: float = np.deg2rad(rec_lon - sp_lon)
    S: float = np.sqrt(np.sum((sp_loc - rec_loc) ** 2))
    Z: float = angle_between_vectors(vec1=[0, 0, rec_loc[2]], vec2=np.array(rec_loc) - np.array(sp_loc))
    return [S * np.sin(Z) * np.cos(A), S * np.sin(Z) * np.sin(A), S * np.cos(Z)]


class era5:
    def __init__(self, path_ws_t: str, path_swh: str):
        era_ws_t: nc.Dataset = nc.Dataset(path_ws_t)
        era_swh: nc.Dataset = nc.Dataset(path_swh)
        self.era_time: list = era_ws_t["valid_time"][:].tolist()
        self.month_day: list = [(i * 24, i * 24 + 24) for i in range(0, int(len(self.era_time) / 24))]
        self.era_ws_lat: list =  era_ws_t["latitude"][:]
        self.era_ws_lon: list = era_ws_t["longitude"][:]
        self.era_swh_lat: list = era_swh["latitude"][:]
        self.era_swh_lon: list = era_swh["longitude"][:]
        self.len_swh_lat:int = len(self.era_swh_lat)
        self.len_swh_lon:int = len(self.era_swh_lon)
        self.len_ws_lat:int = len(self.era_ws_lat)
        self.len_ws_lon:int = len(self.era_ws_lon)
        self.uws: np.ma.core.MaskedArray = era_ws_t["u10"][:]
        self.vws: np.ma.core.MaskedArray = era_ws_t["v10"][:]
        self.skt: np.ma.core.MaskedArray = era_ws_t["skt"][:]
        self.iws: np.ma.core.MaskedArray = era_ws_t["i10fg"][:]
        self.sst: np.ma.core.MaskedArray = era_ws_t["sst"][:]
        self.swh: np.ma.core.MaskedArray = era_swh["swh"][:]
        self.uws100:np.ma.core.MaskedArray = era_ws_t["u100"][:]
        self.vws100:np.ma.core.MaskedArray = era_ws_t["v100"][:]
        del era_swh, era_ws_t


class cyL1data:
    def __init__(self, file_path: str):
        self.name: str = file_path.split("/")[-1]
        self.year: int = int(self.name[12:16])
        self.month: int = int(self.name[16:18])
        self.day: int = int(self.name[18:20])
        self.cynum: int = int(self.name[3:5])
        self.data: nc.Dataset = nc.Dataset(file_path)
        self.nbrcs: np.array = np.array(self.data["ddm_nbrcs"])
        self.les: np.array = np.array(self.data["ddm_les"])
        self.lat: np.array = np.array(self.data["sp_lat"])
        self.lon: np.array = np.array(self.data["sp_lon"])
        self.snr: np.array = np.array(self.data["ddm_snr"])
        self.qc: np.array = np.array(self.data["quality_flags"])
        self.qc2: np.array = np.array(self.data["quality_flags_2"])
        self.sp_rx_gain: np.array = np.array(self.data["sp_rx_gain"])
        self.ddm_brcs_uncert: np.array = np.array(self.data["ddm_brcs_uncert"])
        self.inc: np.array = np.array(self.data["sp_inc_angle"])
        self.time:np.array = np.array(self.data["ddm_timestamp_utc"])
        self.sp_pos_x:np.array = np.array(self.data["sp_pos_x"])
        self.sp_pos_y:np.array = np.array(self.data["sp_pos_y"])
        self.sp_pos_z:np.array = np.array(self.data["sp_pos_z"])
        self.tx_vel_x:np.array = np.array(self.data["tx_vel_x"])
        self.tx_vel_y:np.array = np.array(self.data["tx_vel_y"])
        self.tx_vel_z:np.array = np.array(self.data["tx_vel_z"])
        self.sc_pos_x:np.array = np.array(self.data["sc_pos_x"])
        self.sc_pos_y:np.array = np.array(self.data["sc_pos_y"])
        self.sc_pos_z:np.array=np.array(self.data["sc_pos_z"])
        self.sc_lon:np.array=np.array(self.data["sc_lon"])
        self.sc_vel_x:np.array=np.array(self.data["sc_vel_x"])
        self.sc_vel_y:np.array = np.array(self.data["sc_vel_y"])
        self.sc_vel_z:np.array=np.array(self.data["sc_vel_z"])
        self.prn:np.array = np.array(self.data["prn_code"])
        self.sp_az_body:np.array = np.array(self.data["sp_az_body"])
        self.sp_az_orbit:np.array = np.array(self.data["sp_az_orbit"])
        self.brcs :np.array = np.array(self.data["brcs"])
        self.brcs_ddm_sp_bin_delay_row:np.array = np.array(self.data["brcs_ddm_sp_bin_delay_row"])
        self.brcs_ddm_sp_bin_dopp_col:np.array=np.array(self.data["brcs_ddm_sp_bin_dopp_col"])
        self.eirp = np.array(self.data["gps_eirp"])
        self.power_analog=np.array(self.data["power_analog"])
        self.rx_to_sp_range = np.array(self.data["rx_to_sp_range"])
        self.tx_to_sp_range=np.array(self.data["tx_to_sp_range"])
        self.gps_eirp = np.array(self.data["gps_eirp"])



        del self.data

def get_EdgeA(ddm:np.array or np.ma.core.MaskedArray)->float:
    try:
        ddm = mt.normalization(ddm)
        noise = ddm[:4]
        edge = ddm[4:]
        return np.average(edge)-np.average(noise)
    except:
        return 0

def get_wd_parameter(ddm: np.array, brcs_ddm_sp_bin_delay_row: float, brcs_ddm_sp_bin_dopp_col: float) -> tuple:
    """
    计算与风向相关的几何方位角
    args：
        ddm:power_analog,brcs,raw_counts都可以
        brcs_ddm_sp_bin_delay_row:双基地雷达截面DDM（brcs）镜面点延迟的零基延迟行。请注意，这是一个浮点值
        brcs_ddm_sp_bin_dopp_col:双基地雷达截面DDM（brcs）中镜面点多普勒的零基多普勒列。请注意，这是一个浮点值
    Return：
        两个几何方位角
    """
    try:
        ddm = np.array(ddm)
        ddm_norm: np.array = np.array((ddm - np.min(ddm)) / (np.max(ddm) - np.min(ddm))).T
        mid = np.where(ddm_norm>(1/np.e),1,0)
        ddm_norm = mid*ddm_norm

        ddm_sk: np.array = np.array(
            [[ddm_norm[i][j] if 0.3 < ddm_norm[i][j] < 0.7 else 0 for j in range(0, len(ddm_norm[i]))] for i in
             range(0, len(ddm_norm))])


        ddm_cm: tuple = get_centroid(ddm_norm)
        ddm_sk_cm: tuple = get_centroid(ddm_sk)
        return (az(brcs_ddm_sp_bin_delay_row, brcs_ddm_sp_bin_dopp_col, ddm_cm[1], ddm_cm[0]),
                az(brcs_ddm_sp_bin_delay_row, brcs_ddm_sp_bin_dopp_col, ddm_sk_cm[1], ddm_sk_cm[0]))
    except Exception as e:
        print(f"捕获到异常:{e}")
        raise
def get_pr(ddm: np.array, brcs_ddm_sp_bin_delay_row: float, brcs_ddm_sp_bin_dopp_col: float) -> float:
    """
    计算(power ratio,PR)功率比
    args:
        ddm:power_analog,brcs,raw_counts 都可以
        brcs_ddm_sp_bin_delay_row:双基地雷达截面DDM（brcs）镜面点延迟的零基延迟行。请注意，这是一个浮点值
        brcs_ddm_sp_bin_dopp_col:双基地雷达截面DDM（brcs）中镜面点多普勒的零基多普勒列。请注意，这是一个浮点值
    Return：
        PR
    """
    try:
        delay = round(brcs_ddm_sp_bin_delay_row)
        doppler = round(brcs_ddm_sp_bin_dopp_col)
        sum_ddma = np.sum(ddm[delay:delay + 3, doppler:doppler + 3])
        ddmout = np.sum(ddm) - sum_ddma
        return sum_ddma / ddmout
    except Exception as e:
        print(f"捕获到异常:{e}")
        raise
def get_Moment_of_inertia(ddm: np.array or np.ma.core.MaskedArray) -> float:
    """
    计算归一化转动惯量
    args:
        ddm:power_analog,raw_counts,brcs都可以
    Return:
        二维归一化转动惯量
    """
    try:
        m :float= np.sum(ddm)
        gx :float = np.sum([[x * ddm[x][y] for y in range(0, len(ddm[x]))] for x in range(0, len(ddm))]) / m
        gy:float = np.sum([[y * ddm[x][y] for y in range(0, len(ddm[x]))] for x in range(0, len(ddm))]) / m
        J :float = np.sum(
            [[((x - gx) ** 2 + (y - gy) ** 2) * ddm[x][y] for y in range(0, len(ddm[x]))] for x in range(0, len(ddm))])
        return np.sqrt(J) / m
    except Exception as e:
        print(f"捕获到异常:{e}")
        raise
def get_sr(power_analog: np.array, rxtosp: float, txtosp: float, gps_eirp: float, sp_rx_gain: float,
           inc: float) -> float:
    """
    计算地表反射率SR
    args:
        power_analog:DDM仓功率值，17×11
        rxtosp:rx_to_sp_range,CYGNSS卫星与镜面反射点之间的距离
        txtosp:tx_to_sp_range,GPS卫星与镜面反射点之间的距离
        gps_eirp:GPS有效各向同性辐射功率
        sp_rx_gain:镜面反射点方向上的接收天线增益
        inc:sp_inc_angle,信号入射角
    Return:
        地表反射率SR
    """
    try:
        pcrl: float = np.max(power_analog)
        s_r = (np.log10(pcrl) + 2 * np.log10(4 * np.pi) +
               2 * np.log10(rxtosp + txtosp) -
               np.log10(gps_eirp) - np.log(sp_rx_gain) - 2 * np.log(0.19)) * 10
        sr = s_r - 10 * np.log10(np.cos(inc * (np.pi / 180)))
        return sr
    except Exception as e:
        print(f"捕获到异常:{e}")
        raise
def get_centroid(ddm: np.array) -> list or tuple:
    """
    计算质心
    args:
        ddm:power_analog,raw_counts,brcs都可以
    Return:
        二维质心坐标
    """
    try:
        return (np.sum([[x * ddm[x][y] for y in range(0, len(ddm[x]))] for x in range(0, len(ddm))]) / np.sum(ddm),
                np.sum([[y * ddm[x][y] for y in range(0, len(ddm[x]))] for x in range(0, len(ddm))]) / np.sum(ddm))
    except Exception as e:
        print(f"捕获到异常:{e}")
        raise

def convert_time(time: int) -> datetime.datetime:
    return datetime.datetime(year=1970, month=1, day=1, hour=0, minute=0) + datetime.timedelta(seconds=time)


def match_data(cy_file_path: str, era_ws_path: str, era_swh_path: str,output:str ):
    era_data: era5 = era5(path_ws_t=era_ws_path, path_swh=era_swh_path)

    cy_day_list: list = os.listdir(cy_file_path)
    for constellation in cy_day_list[5:]:
        """
        constellation 为日期中的日
        """

        satellites: list = os.listdir("/".join([cy_file_path, constellation]))
        day: int = int(satellites[0][18:20]) - 1
        miduws: np.ma.core.MaskedArray = era_data.uws[
                                         era_data.month_day[day][0]:era_data.month_day[day][1]]
        midvws: np.ma.core.MaskedArray = era_data.vws[
                                         era_data.month_day[day][0]:era_data.month_day[day][1]]
        midiws: np.ma.core.MaskedArray = era_data.iws[
                                         era_data.month_day[day][0]:era_data.month_day[day][1]]
        midsst: np.ma.core.MaskedArray = era_data.sst[
                                         era_data.month_day[day][0]:era_data.month_day[day][1]]
        midskt: np.ma.core.MaskedArray = era_data.skt[
                                         era_data.month_day[day][0]:era_data.month_day[day][1]]
        midswh: np.ma.core.MaskedArray = era_data.swh[
                                         era_data.month_day[day][0]:era_data.month_day[day][1]]
        miduws100: np.ma.core.MaskedArray = era_data.uws100[
                                         era_data.month_day[day][0]:era_data.month_day[day][1]]
        midvws100: np.ma.core.MaskedArray = era_data.vws100[
                                            era_data.month_day[day][0]:era_data.month_day[day][1]]
        out:matched = matched()
        a=time.time()
        for satellite in satellites:
            cydata: cyL1data = cyL1data(file_path="/".join([cy_file_path, constellation, satellite]))
            c1=np.where(cydata.les > 0, 0, 1)
            c2=np.where(cydata.nbrcs > 0, 0, 1)
            c3=np.where(cydata.qc == 0, 0, 1)
            c4=np.where(cydata.qc2 == 0, 0, 1)
            c5=np.where(cydata.inc < 60, 0, 1)
            c6=np.where(cydata.lat > -40, 0, 1)
            c7=np.where(cydata.lat < 40, 0, 1)
            control_index: np.array = c1+c2+c4+c5+c6+c7

            for sample in range(len(control_index)):
                for channel in range(0, 4):
                    if control_index[sample][channel]:
                        pass
                    else:
                        """
                        aa = round(a) if a>=0 else round(360+a)
                        """
                        cy_sp_lat: float = cydata.lat[sample][channel]
                        cy_sp_lon =cydata.lon[sample][channel]
                        cy_sp_lon_c: float = cy_sp_lon if cy_sp_lon<=180 else cy_sp_lon-360
                        swh_lat_index: int = era_data.len_swh_lat-int(round((cy_sp_lat +40)/0.5))-1
                        swh_lon_index: int = int(round(cy_sp_lon_c/0.5))+180*2
                        swh_lon_index:int = swh_lon_index if swh_lon_index<era_data.len_swh_lon else 0

                        ws_lat_index :int = era_data.len_ws_lat-int(round((cy_sp_lat+40)/0.25))-1
                        ws_lon_index:int = int(round(cy_sp_lon_c/0.25))+180*4
                        ws_lon_index:int = ws_lon_index if ws_lon_index<era_data.len_swh_lon else 0

                        time_index_mid:int = int(np.round(cydata.time[sample]/3600))
                        time_index :int = time_index_mid if time_index_mid<24 else 23
                        uws:float = miduws[time_index][ws_lat_index][ws_lon_index]
                        vws:float = midvws[time_index][ws_lat_index][ws_lon_index]
                        uws100: float = miduws100[time_index][ws_lat_index][ws_lon_index]
                        vws100: float = midvws100[time_index][ws_lat_index][ws_lon_index]
                        mid_loc = topocentric_coordinate_system(sp_loc=np.array([cydata.sp_pos_x[sample][channel],
                                                                                 cydata.sp_pos_y[sample][channel],
                                                                                 cydata.sp_pos_z[sample][channel], ], ),
                                                                rec_loc=np.array([cydata.sc_pos_x[sample],
                                                                                  cydata.sc_pos_y[sample],
                                                                                  cydata.sc_pos_z[sample]]),
                                                                sp_lon=cy_sp_lon,
                                                                rec_lon=cydata.sc_lon[sample])
                        sr = get_sr(power_analog=cydata.power_analog[sample][channel],
                                    rxtosp=cydata.rx_to_sp_range[sample][channel],
                                    txtosp=cydata.tx_to_sp_range[sample][channel],
                                    gps_eirp=cydata.gps_eirp[sample][channel],
                                    sp_rx_gain=cydata.sp_rx_gain[sample][channel],
                                    inc=cydata.inc[sample][channel]
                                    )
                        pr = get_pr(ddm=cydata.brcs[sample][channel],
                                    brcs_ddm_sp_bin_delay_row=cydata.brcs_ddm_sp_bin_delay_row[sample][channel],
                                    brcs_ddm_sp_bin_dopp_col=cydata.brcs_ddm_sp_bin_dopp_col[sample][channel])
                        nmi = get_Moment_of_inertia(ddm=cydata.brcs[sample][channel])
                        aziele: list = spherical_coordinates(x=mid_loc[0], y=mid_loc[1], z=mid_loc[2])
                        out.content["azi"].append(aziele[0])
                        out.content["ele"].append(aziele[1])
                        out.content["ews"].append(np.sqrt(uws**2+vws**2))
                        out.content["ewsu"].append(uws)
                        out.content["ewsv"].append(vws)
                        out.content["ews100"].append(np.sqrt(uws100 ** 2 + vws100 ** 2))
                        out.content["ews100u"].append(uws100)
                        out.content["ews100v"].append(vws100)
                        out.content["iws"].append(midiws[time_index][ws_lat_index][ws_lon_index])
                        out.content["sst"].append(midsst[time_index][ws_lat_index][ws_lon_index])
                        out.content["skt"].append(midskt[time_index][ws_lat_index][ws_lon_index])
                        out.content["swh"].append(midswh[time_index][swh_lat_index][swh_lon_index])
                        out.content["lat"].append(cy_sp_lat)
                        out.content["lon"].append(cy_sp_lon)
                        out.content["nbrcs"].append(cydata.nbrcs[sample][channel])
                        out.content["les"].append(cydata.les[sample][channel])
                        out.content["snr"].append(cydata.snr[sample][channel])
                        out.content["inc"].append(cydata.inc[sample][channel])
                        out.content["cynum"].append(cydata.cynum)
                        out.content["channel"].append(channel)
                        out.content["txvx"].append(cydata.tx_vel_x[sample][channel])
                        out.content["txvy"].append(cydata.tx_vel_y[sample][channel])
                        out.content["txvz"].append(cydata.tx_vel_z[sample][channel])
                        out.content["scvx"].append(cydata.sc_vel_x[sample])
                        out.content["scvy"].append(cydata.sc_vel_y[sample])
                        out.content["scvz"].append(cydata.sc_vel_z[sample])
                        out.content["prn"].append(cydata.prn[sample][channel])
                        out.content["sp_rx_gain"].append(cydata.sp_rx_gain[sample][channel])
                        out.content["sp_az_body"].append(cydata.sp_az_body[sample][channel])
                        out.content["sp_az_orbit"].append(cydata.sp_az_orbit[sample][channel])
                        
                        azlist:tuple = get_wd_parameter(ddm=cydata.brcs[sample][channel],
                                                        brcs_ddm_sp_bin_delay_row=cydata.brcs_ddm_sp_bin_delay_row[sample][channel],
                                                        brcs_ddm_sp_bin_dopp_col=cydata.brcs_ddm_sp_bin_dopp_col[sample][channel])
                        out.content["a1"].append(azlist[0])
                        out.content["a2"].append(azlist[1])
                        out.content["edgea"].append(get_EdgeA(cydata.brcs[sample][channel]))
                        out.content["sr"].append(sr)
                        out.content["pr"].append(pr)
                        out.content["nmi"].append(nmi)
                        out.content["eirp"].append(cydata.eirp[sample][channel])
                        out.content["time"].append(cydata.time[sample])

        out.save(output+cy_file_path.split("/")[-2]+constellation)
        b=time.time()
        print(output+constellation,b-a)






match_data(cy_file_path="G:/CYdata/L1V3.1/2023/02/",
           era_ws_path="G:/ERA-5_new/2023/ERA-5_202302_ws_t.nc",
           era_swh_path="G:/ERA-5_new/2023/ERA-5_202302_swh.nc",
           output="G:/nwnw/2023/02/")



